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Abstract. Magma viscosity is strongly temperature-dependent. 
When hot magma flows in a conduit, heat is lost through the walls 
and the temperature decreases along the flow causing a viscosity 
increase. For particular values of the controlling parameters the 
steady-flow regime in a conduit shows two stable solutions belong- 
ing either to the slow or to the fast branch. As a consequence, this 
system may show an hysteresis effect, and the transition between 
the two branches can occur quickly when certain critical points are 
reached. In this paper we describe a model to study the relation be- 
tween the pressure at the inlet and the volumetric magma flow rate in 
a conduit. We apply this model to explain an hysteric jump observed 
during the dome growth at Soufriere Hills volcano (Montserrat), and 
described by Melnik and Sparks [1999] using a different model. 



1. Introduction 

Like other liquids of practical interest (e.g. polymers, oils), 
magma has a strong temperature-dependent viscosity. It is well 
known that the flow in conduits of this kind of fluids admits one or 
more solutions for different ranges of the controlling parameters. 
The problem of the multiple solutions and their stability was stud- 
ied by Pearson et al. [1973] and by Aleksanopol'skii and Naidenov 
[1979]. More recently this phenomenon, applied to magma flows, 
was investigated with more sophisticated models by Helfrich [1995] 
and by Wylie and Lister [1995]. In Skid'skiy et al. [1999] a simple 
method to find the conditions of multiple solutions was adopted and 
an experimental study to verify the hysteresis phenomenon predicted 
for this process was performed. 

2. The model 

In this paper we investigate a simple one-dimensional flow model 
of a fluid with temperature-dependent viscosity, with the essential 
physical properties characterizing the phenomenon of the multiple 
solutions and hysteresis as m Skid'skiy etal. [1999]. The fluid flows 
in a conduit with constant cross section and constant temperature at 
the wall boundaries. We assume that the fluid properties are constant 
with the temperature except for the viscosity, and we neglect the heat 
conduction along the streamlines, and the viscous heat generation. 
Moreover, we assume a linear relation between the shear stress and 



^AIso at Dipartimento di Scienze della Terra, Universita degli Studi di 
Pisa, Italy. 



Copyright 2002 by the American Geophysical Union. 



Paper number 2001GL014493. 
0094-8276/08/2001GL014493$5.00 



the strain rate (Newtonian rheology). This last assumption is intro- 
duced to simplify the model and allows us to demonstrate that the 
multiple solutions in conduit flows are the direct consequence of 
the viscosity increase along the conduit induced by cooling (under 
particular boundary conditions). Under these hypotheses, the equa- 
tions for momentum and energy balance for the one-dimensional 
steady flow in a long circular conduit (R <C L) at low Reynolds 
number are: 



1_9 f^^j,yd}i\ ^ dP 
r dr \ dr J dz 

(kr — \ - c V — 
r dr \ dr J ^ ^ dz 



(1) 
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where R is the conduit radius, L conduit length, r radial direction, 
z direction along the flow, v velocity along the flow (we assume 
V — v{r, z) ~ v{r)), Cp specific heat at constant pressure, k ther- 
mal conductivity, p density, /i viscosity, and P pressure, and T 
temperature. For magma, the dependence of the viscosity on the 
temperature is well described by the Arrhenius law: 



fi ^ fj,A exp{B/T) 



(3) 



where fiA is a constant and B the activation energy. In this paper, 
for simplicity, we approximate eq. (3) by the Nahme's exponential 
law, valid when (T — Tr)/Tr <C 1, where Tr is a reference tem- 
perature: 



fj,^ fj,R exp [~(3{T — Tr)] 



(4) 



with /3 — B /Tr and /ir — fiA exp{B /Tr). Following Skid'skiy 
et al. [1999], we introduce two new variables: the volumetric flow 
rate 



Q = 2iT v{r)rdr 



and the convected mean temperature 

2ir 



T*{z) 



T{r, z)v{r)rdr 



(5) 



(6) 



To satisfy the mass conservation, the volumetric flow rate Q is con- 
stant along the flow. Integrating eq. (1) and (2) and expressing the 
solutions in terms of Q and T*, we obtain: 



ttP" exp [I3{T* - T^ 

^dT* 

pCpQ— 



dP 



dz 



2-KRa{T^ - T*) 



(7) 



(8) 



where T^ is the wall temperature (taken here as reference, T^ = 
Tr), (Xw fluid viscosity at T = T^, and a the coefficient of heat 
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transfer through the walls defined as: 

Tu,-T* dr , 
In the averaged model (eq. 7), we have adopted 

M«/iiiexp[-/3(T* -Tfl)] 



(9) 



(10) 



and as pointed out in Helfrich [1995], the value of /3 in eq. (10), is 
usually smaller than the actual value of /3 in the fluid (eq. 4). 
Eq. (7) and (8) with the boundary conditions: 



P(0) = Po, P{L)=Q, T*(0) = To 



(11) 



give an approximate solution of the problem. These equations are 
similar to the model of Pearson et al. [1973] (for a plane flow) and 
were used by Shah and Pearson [1974] to study the viscous heating 
effects. In the nondimensional form, eq. (7) and (8) are: 



dp 

dC 

dd 



= -qe 



(12) 



where 



9 = 



2nRLa ' 

pCpR^ 



(13) 



The boundary conditions (1 1) are rewritten for the new variables: 



p{Q) =po = 



p(l) = 0, e{Q)=B (14) 



withB = i3{To - r»). The solutions of eq. (12), satisfying eq. (14) 
at the boundaries, are: 




Figure 1 

Figure 1. Plot of the relation between the nondimensional flow rate 
q and the nondimensional pressure at conduit inlet p, for different 
values of the narameter K. resulting from ea. n6V 
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Figure 2. Plot of the experimental results (crosses) of Skul 'skiy etal. 
[1999]. The dashed line follows the "history" of the flow regimes 
imposed to the fluid; the full line represents the path predicted by 
the model with the same parameters of the experiment. Modified 
aStet Skul'skiy et al. [1999]. 



P(C)-P0 = g/°exp(-5e-f/''K 
e(C) = Bexp(-C/g) 



(15) 



and, therefore, the relation between the nondimensional pressure 
at the conduit inlet po and the nondimensional flow rate q is: 

po^q [ exp(-i3e-'^/«)dC (16) 
Jo 

In Fig. 1 we plot relation (16) obtained numerically. We observe 
that for values of B greater than a critical value Be ^ 3, there are 
values of po which correspond to three different values of q. By us- 
ing a simpler model, Skid 'skiy et al. [1999] found Be = 4, whereas 
Helfrich [1995] found Be = 3.03, in good agreement with Pear- 
son et al. [1973]. Moreover, Shah and Pearson [1974] showed that 
when the viscous heat generation is important, the value of Be can 
be lower, but for high values of B, the relation between po and q is 
similar to the case without viscous heat generation. 

The stabiUty analysis of the three branches (slow, intermediate 
and fast) shows that the intermediate branch is never stable. More- 
over, one part of the slow branch is stable to two-dimensional pertur- 
bations but unstable to the three-dimensional ones, in a way similar 
to the Saff man-Taylor instability [Wylie and Lister, 1995]. In the 
case of multiple solutions, an hysteresis phenomenon may occur, as 
proposed in Wylie and Lister [1995] and verified experimentally in 
Skul'skiy et al. [1999]. 

In the experiments of Skul'skiy et al. [1999] a fluid with pre- 
scribed temperature and pressure was injected into a capillary tube 



Table 1. Typical Values of the Parameters Used in This Paper 



Parameter 


Symbol 


Value 


Unit 


Rock density 


Pr 


2600 


Kgm-3 


Magma density 


Pm 


2300 


Kgm-3 


Conduit radius 


R 


15 


m 


Conduit length 


L 


5000 


m 


Magma temperature 


To 


1123 


K 


Magma specific heat 


Cp 


1000 


JKg-i K-i 


Magma thermal conductivity 


k 


2 


Wm-i K-i 
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with constant wall temperature and controlled fluid pressure and 
flow rate. The device is used to show the transition between the two 
regimes corresponding to the upper and the lower branches. 

A comparison between the experimental results (crosses) and the- 
ory (full line) is shown in Fig. 2 for the nondimensional variables po 
and q, for B = 4.6. The dashed lines indicate the pressure history 
prescribed in the experiments. The two steady-state regimes corre- 
sponding to the slow and to the fast branch were clearly recorded. 
Starting with a low pressure conflguration (point A) and by increas- 
ing the pressure, the flow rate increases along the slow branch until 
it reaches a critical point (point B). Here, a jump to the fast branch 
occurs (point C). Increasing the pressure further, the flow rate in- 
creases along this branch, whereas, by decreasing the pressure, the 
flow rate decreases moving along the upper branch, until it reaches 
another critical point (point D) where the jump on the slow branch 
occurs (point E). On the slow branch the nondimensional flow rate is 
more than one order of magnitude lower than that on the fast branch. 



3. Application 



A t5T5ical value of the viscosity is /ito = n{To) = lO'^Pas for 
Montserrat andesite at 1123 K with 4% water content [Melnik and 
Sparks, 1999]. This value is in good agreement with eq. (17); 
in fact, for B = B'"""* and assuming w 873 K, we have 
HO = /it"'*exp(B) ~ 1.6 X lO'^Pas. Moreover, for B = B*""* 
and To - r„ « 250 K gives /3 ~ 0.014 K"^ (for example, from 
data of Hess and Dingwell [ 1 996] for a magma with a similar com- 
position and 4% water content, we obtain /3 ~ 0.016 K~^). Finally, 
for the heat transfer coefficient we have: 

Q;Ri^«4 Wm-^K-i (18) 

Ot 

where 5t is the thermal boundary layer of the flow, while using 
fe = 2 W m"^ K"\ (5t « 50 cm {Bruce and Happen [1989] used 
(5t ~ 10 cm for dyke length of about one kilometer). 

Moreover, we verify the basic assumptions of the model: the 
assumption of one-dimensional flow is based on the small diame- 
ter/lengthratiooftheconduit(J?/I/ ~ 10~^), andthesmallReynold 



During some basaltic fissure eruptions in Hawaii and in Iceland, 
the eruption begins with a rapid opening at high flow rate and, after 
few hours, the flow rate quickly decreases. To explain this phe- 
nomenon, Wylie and Lister [1995] and Helfrich [1995] proposed 
a model similar to the model presented in this paper, based on the 
hysteric jump between the fast branch and the slow branch when 
the pressure driving the eruption decreases. 

A similar phenomenon, showing a jump in the mass flow rate, 
was observed during the dome growth (1995- 1999) at Soufriere Hills 
volcano in Montserrat as described by Melnik and Sparks [1999]. 
The model used by Melnik and Sparks [1999] to explain this ef- 
fect is essentially based on the crystal growth kinetics which affects 
magma viscosity, and the mechanical coupling between the gas and 
the melt through the Darcy law. 

In this paper we explain the same phenomenon in terms of the 
viscosity variation governed by cooling along the flow. However, 
since the crystal content is physically related to the magma temper- 
ature, the two models are physically related. 

Using the data of Tab. 1, we fit the curve of eq. (16) with the 
observed values of the discharge rates and dome height reported in 
Melnik and Sparks [1999]. 

The variation of the dome height reflects a change of the exit 
pressure and, as a consequence, a variation of the difference be- 
tween the inlet and the outlet pressures. Since we assume a zero 
exit pressure in our model and the gravity term was not explicitly 
accounted for in eq. (1), the variable po represents the overpressure 
at the base of the conduit (to obtain the actual value, the hydrostatic 
pressure pg[L + H) must be added to its value, where L is the 
conduit length and H is the dome height). 

The values of a, B and Hw are chosen by least square best fit- 
ting of the observed data, and the wall temperature was defined 
as the temperature for which magma ceases to flow. Fig. 3 shows 
the results of least square fitting: the discharge rate is reported on 
the X-axis and the overpressure on the y-axis. The crosses indicate 
the observed values, reported in Melnik and Sparks [1999], and the 
dome height expressed in terms of overpressure. 

The values obtained by the best fit of the observed data are: 
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Figure 3. Relation between the discharge rate and the pressure. Full 
line represents the model prediction; crosses represents observed 
values at Soufriere Hills from Melnik and Sparks [1999]. 




B 



6est 



B" 



3.56 Wm^^K-i 
5 X 10* Pa s 
3.46 



(17) 



Fig. 3 shows the agreement of the model with the observed data; the 
proposed model is able to explain the hysteresis effect observed in 
dome growth at Soufriere Hill (Montserrat) by Melnik and Sparks 
[1999], and modeled in a different way. 



Figure 4 



Fioiire 4. Possible nulsatin? behavior nredicted hv the model. 
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number is simply verified: 

pRv _ 2300 X 15 X 0.003 5 .^o^ 

-n-e = ~ T-^ 10 (19) 

Mo 10^ 

The viscous heating effects can be neglected because the Nahme 
number based on the shear stress is small: 

g = -^—^ — « 1 (20) 

The assumption of negligible heat conduction along the stream- 
lines is justified by the high value of the Peclet number (the ra- 
tio between the advective and the conductive heat conduction): 

Pe = {pCpVR)/k « 10^ 

The existence of the multiple solutions for the steady flow al- 
lows the system to show a pulsating behavior between the different 
solutions. In the case where the initial pressure, is greater than the 
critical pressure corresponding to point E in Fig. 4, the system is on 
the fast branch of the solution, such as in point A. If the pressure de- 
creases, the system moves along this branch up to the critical point 
C. In C a jump to the slow branch occurs (point D). If the pressure 
continues to decrease, the discharge rate tends to zero. Instead, if 
the pressure increases, the system moves along the slow branch up 
to the other critical point E. In this point the jump occurs on the fast 
branch and the system reaches point B. 

The overpressure conditions and pulsating activity, typical of 
dome eruptions, are evident not only at Soufriere Hills, but also 
in Santiaguito (Guatemala), Mount Unzen (Japan), Lascar (Chile), 
Galeras (Colombia) and Mount St. Helens (USA) [Melnik and 
Sparks, 1999]. 

4. Conclusion 

Magma flow in conduits shows the existence of multiple solu- 
tions, like other fluids with sttong temperature dependent viscosity. 
This is a consequence of the increase of viscosity along the flow 
due to cooling. For a given pressure drop along the conduit, one or 
two stable regimes (fast and slow branches) may exist. The tran- 
sition between the two branches occurs when critical values are 
reached, and an hysteresis phenomenon is possible. These jumps 
were evident during the dome growth in the 1995-1999 Soufriere 
Hills (Montserrat) eruption. The pulsating behaviour of the dome 
growth was previously modeled by Melnik and Sparks [1999] in 



terms of the nonlinear effects of crystalhzation and gas loss by per- 
meable magma. 

In this paper we propose a model to describe the nonlinear jumps 
between the two stable solutions as a consequence of the coupling 
between the momentum and energy equation induced by the strong 
temperature-dependent viscosity of magma. 

However, since the crystal content is a consequence of cooling, 
the two models, although different, are physically related. 
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